0 


PL-TR-97-2150 

Advances  in  Numerical  Methods 


K.  L.  McLaughlin 
B.  Shkoller 
D.  E.  Wilkins 
T.  G.  Barker 


Maxwell  Technologies 
8888  Balboa  Avenue 
San  Diego,  CA  92123-1506 


September1997  19980421  120 

Final  Report 

in  .  n  .  DTIC  QUALITY  INSPECTED  3 

10  August  1995  -  9  August  1997  * 


Approved  for  Public  Release;  Distribution  Unlimited 


U.S.  DEPARTMENT  OF  ENERGY 

Office  of  Nonproliferation  and  National  Security 

WASHINGTON,  DC  20585 

AIR  FORCE  RESEARCH  LABORATORY 
Space  Vehicles  Directorate 
29  Randolph  Road 

AIR  FORCE  MATERIEL  COMMAND 
HANSCOM  AFB,  MA  0173103010 


SPONSORED  BY 
Department  of  Energy 

Office  of  Non-Proliferation  and  National  Security 

MONITORED  BY 
Air  Force  Research  Laboratory 
CONTRACT  No.  F19628-95-C-0207 


The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors  and  should  not  be 
interpreted  as  representing  the  official  policies,  either  express  or  implied,  of  the  Air  Force  or  U.S. 
Government. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


KATHARINE  KADINSKY-CADE 
Contract  Manager 


CHARLES  P.  PIKE,  Deputy  Director 
Integration  and  Operations  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is  releasable  to  the  National 
Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  copies  from  the  Defense  Technical  Information  Center.  All  others  should 
apply  to  the  National  Technical  Information  Service. 


If  your  address  has  changed,  or  you  wish  to  be  removed  from  the  mailing  list,  or  if  the  addressee  is  no 
longer  employed  by  your  organization,  please  notify  AFRL/VSOE,  29  Randolph  Road,  Hanscom  AFB, 
MA  01731-3010.  This  will  assist  us  in  maintaining  a  current  mailing  list. 


Do  not  return  copies  of  the  report  unless  contractual  obligations  or  notices  on  a  specific  document  requires 
that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  ia  estimated  to  average  1  hour  par  response,  including  the  time  for  reviewing  instruction*,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704-0188),  Washington.  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

September  1 997  Final  -  08/10/95-08/09/97 

4.  TITLE  AND  SUBTITLE 

Advances  in  Numerical  Methods 

5.  FUNDING  NUMBERS 

PE69120H 

PR  DENNTAGM 

WU  AT 

FI  9628-95-C-0207 

6.  AUTHOR(S) 

Keith  L.  McLaughlin  D.E.  Wilkins 

B.  Shkoller  T.G.  Barker 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Maxwell  Technologies 

8888  Balboa  Avenue 

San  Diego,  CA  92123-1506 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

MFD-TR-97-15870 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Road 

Hanscom  AFB,  MA  01 731  -301 0 

Contract  Manager:  Katherine  Kadinsky-Cade/VSBS 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

PL-TR-97-2150 

11.  SUPPLEMENTARY  NOTES 

This  research  was  sponsored  by  the  Department  of  Energy,  Office  of  Non-Proliferation  and 
National  Security,  Washington,  DC  20585. 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

This  report  summarizes  work  conducted  over  a  two  year  period  to  improve  numerical 
methods.  The  report  has  three  main  sections.  Section  2  covers  work  on  synthetics  for 
layered  earth  models.  Section  3  reports  on  three-dimensional  (3D)  elastic  finite  difference 
calculations  and  the  relationship  between  one-dimensional  (ID)  (layered  earth)  and  3D 
laterally  heterogeneous  wavefields.  Section  4  reports  on  development  of  a  finite  difference 
program  to  model  the  response  of  permeable  hoses  to  atmospheric  pressure  fluctuations 
from  wind  turbulence. 

14.  SUBJECT  TERMS 

Synthetic  Seismograms  3D  Finite  Differences 

Infrasound  Permeable  Hoses 

15.  NUMBER^  PAGES 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 

20.  LIMITATION  OF  ABSTRACT 

SAR 

NSN  7540-01-280-5500  Standard  Form  298  {Rev.  2-89) 


Proscribed  by  ANSI  Std.  Z39-18 


298-102 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


n 


Table  of  Contents 


1.0  Introduction . 1 

2.0  Layered  Earth  Model  Wavenumber  Integration  Synthetics . 2 

2. 1  Wavenumber  Integration  Synthetics  Using  A  Parallel  Virtual  Machine . . 2 

2.1.1  What  is  PVM? . 2 

2.1.2  Wavenumber  Integration  Synthesis  -“The  Perfect  Parallel  Algorithm” . 3 

2.2  Random  Layering . 5 

2.3  Frequency  Dependent  Attenuation,  Q(f)  =  Qo*^ . 11 

2.4  Acknowledgments . 16 

2.5  References . , . 16 

3.0  Finite  Difference  3D  Regional  Scattering  Calculations . 17 

3.1  Introduction . 17 

3.2  Modal  Spectra . 24 

3.3  Wavenumber  Spectra..., . 29 

3.4  Conclusions . . . 33 

3.5  Acknowledgements . 34 

3.6  References . 34 

4.0  Permeable  Hose  Characteristics  and  Noise  Reduction  For  Infrasound 

Monitoring . 35 

4.1  Summary . 35 

4.2  Introduction . 35 

4.3  Theory  of  Noise  Reduction  Hoses  and  Their  Calibration . 37 

4.4  Analysis  of  a  Helmholtz  Oscillator . 43 

4.5  Finite  Difference  Modeling,  Maxhose . 44 

4.6  Stability  Conditions . 48 

4.7  Hose  Calibration  Simulation . 49 

4.8  Finite  Difference  Simulations  of  Wind  Noise  Response . 49 

4.9  Conclusions . 52 

4. 10  Acknowledgments . 53 

4. 1 1  References . 53 


iii 


List  of  Figures 


1.  The  master  program  sends  and  receives  messages  through  the  PVM  daemon 
(PVMD)  which  routes  the  messages  to  other  daemons  on  the  way  to  and  from 


slave  modules . 4 

2.  Comparison  of  synthetic  transverse  component  Green's  function,  Gyxy,  at  400 

km,  source  depth  h=15  km,  with  (top)  and  without  (below)  5%  random 
velocity  variations  in  the  layered  model  D2 . 6 

3.  Comparison  of  vertical  component  Green's  function,  G*,  at  400  km,  source 
depth  h=T  km,  with  (top)  and  without  (below)  5%  random  velocity  variations 

in  the  layered  model  D2 . 6 

4.  P-wave  velocities  versus  depth  for  layered  models  DA,  GA,  H2,  and  N4  with 

0%,  5%,  7.5%,  and  10%  velocity  variations . . 7 

5.  Comparison  of  Green's  functions  for  model  DA  with  0%,  2.5%,  5.0%,  7.5%, 

and  1 0%  velocity  variation  (top  to  bottom) . 8 

6.  Comparison  of  Green's  functions  for  model  GA  with  0%,  2.5%,  5.0%,  7.5%, 

and  10%  velocity  variation  (top  to  bottom) . . 8 

7.  Comparison  of  Green's  functions  for  model  H2  with  0%,  2.5%,  5.0%,  7.5%, 

and  10%  velocity  variation  (top  to  bottom) . . . 9 

8.  Comparison  of  Green's  functions  for  model  N4  with  0%,  2.5%,  5.0%,  7.5%, 

and  10%  velocity  variation  (top  to  bottom) . . . 9 

9.  Vertical  seismogram  at  600  km  from  an  earthquake  source  (h=12.5  km)  for 

layered  model  DA  with  0%,  2.5%,  5.0%,  7.5%,  and  10%  velocity  variation . 10 

10.  Vertical  component  Lg/Pg  ratio  at  600  km  from  an  earthquake  source 
(h=12.5  km)  for  layered  model  DA  with  0%,  2.5%,  5.0%,  7.5%,  and  10% 
velocity  variation. . . . 10 


11.  Vertical  component  Lg/Pg  spectral  ratios  for  models  DA,  GA,  H2,  and  N4  for 

an  explosion  source  (h=l  km)  and  0,  2.5,  5,  7.5,  and  10%  RMS  random 
velocity  layers . . . . . 

12.  Compilation  of  average  explosion  Lg/Pg  spectral  ratios  from  Bennett  et  al. 


(1997)  for  explosions  in  Asia  and  Western  US . 13 

13.  Compilation  of  average  Lg/Pg  spectral  ratios  from  Bennett  et  al.  (1997)  for 

earthquakes  in  Asia,  Western  US,  and  Eastern  US . 13 

14.  Comparison  of  vertical  component  synthetics  (high  pass  at  0.5  Hz)  at  200  km 

for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle),  Q3 
(bottom) . . . . . 14 

* 


IV 


15.  Comparison  of  vertical  component  synthetics  (high  pass  at  0.5  Hz)  at  200  km 

for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle),  Q3 
(bottom) . 14 

16.  Comparison  of  vertical  component  synthetics  (high  pass  at  3.0  Hz)  at  200  km 

for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle),  Q3 
(bottom) . 15 

17.  Comparison  of  vertical  component  synthetics  (high  pass  at  3.0  Hz)  at  200  km 

for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle),  Q3 
(bottom) . 15 

18.  Diagram  of  3  levels  of  nested  grids;  1  grid  on  level  1,  4  grids  on  level  2,  and 

16  grids  on  level  3  for  a  total  of  21  grids.  Each  grid  contains  64  x  64  x  128  = 
524,288  cells  for  a  total  of  1 1  million  cells . 19 

19.  Plot  of  the  test  structure  used  in  all  computations.  P-  and  S-wave  velocities  in 

the  crust  are  simple  gradients  with  Poisson's  ratio  of  0.25 . 20 

20.  Shear  velocity  along  a  Y  =  Constant  slice  of  the  model . 20 

21.  Visualization  of  the  3D  random  variation.  Isosurfaces  with  absolute  variation 

greater  ±la  are  shaded  in  the  top  image.  Regions  within  ±la  are  transparent . 21 

22.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  5 

seconds . . 22 

23.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  10 

seconds . 22 

24.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  15 

seconds . 22 

25.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  20 

seconds  . 23 

26.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  25 

seconds . 23 

27.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  30 

seconds . 23 

28.  Comparison  of  finite  difference  and  wavenumber  integration  synthetics  for 
Run  #1  (5%  RMS,  explosion  source)  at  100  km  for  lowpass  filters  0.2,  0.4, 

0.6,  and  0.8  Hz  (top  to  bottom) . 24 

29.  Plots  of  P-SV  modes  at  0.5  Hz  (right)  and  0.75  Hz  (left)  for  the  test  structure 

shown  in  Figure  19 . 26 

30.  Diagram  of  Modal  Spectra.  C(J,F).  C(J,F)  =  0  to  the  right  of  the  modal  cut- 

offline  . 26 

3 1 .  Modal  spectra  at  120  km  (top)  and  60  km  (bottom)  for  models  with  0%  RMS 

(left)  and  5%  RMS  (right)  lateral  heterogeneity . 27 


V 


32.  Modal  Spectral  Ratio  at  120  km  5%  RMS  model  /  0%  RMS  model.  Yellow 

and  Red  shades  are  modes  enhanced  w.r.t.  the  background  model . 28 

33.  Ratio  of  120  km  to  60  km . . . 28 

34.  Diagram  of  wavenumber-wavenumber  spectra  at  frequency,  f. . 30 

35.  Phase  plane  wavenumber  impulse  response  for  a  cylindrical  wave,  120  km 

from  the  source  at  0.4  (left)  and  0.6  Hz  (right) . 30 

36.  Phase  plane  wavenumber  impulse  response  for  an  off-axis  cylindrical  wave, 

.  120  km  from  the  source  at  0.6  Hz . . . 31 

37.  Wavenumber  spectra  of  the  vertical  component  velocity  field  on  a  Y  =  120 

km  phase  plane  at  0.6  Hz  for  the  5%  RMS  (right)  and  0%  RMS  (left) . 32 

38.  Wavenumber  spectra  of  the  radial  component  velocity  field  on  a  Y  =  120  km 

phase  plane  at  0.6  Hz  for  the  5%  RMS  (right)  and  0%  RMS  (left) . 32 

39.  Wavenumber  spectra  of  the  transverse  component  velocity  field  on  a  Y  =  120 

km  phase  plane  at  0.6  Hz  for  the  5%  RMS  (right)  and  0%  RMS  (left) . 33 

40.  Contours  of  detection  threshold  in  Log(Yield  in  Kt)  at  90%  probability  for  2 

infrasound  detections  at  stations  of  the  proposed  IMS . 36 

41.  Contours  of  detection  threshold  in  Log  (Yield  in  Kt)  at  90%  probability  for  3 

infrasound  detections  at  stations  of  the  proposed . 37 

42.  Wind  generated  eddies  create  pressure  fluctuations  and  hence  noise . 37 

43 .  A  hose  element:  length  dx,  radius  a,  flow  in  the  x  direction  qx,  flow  in/out  qa . 38 

44.  A  simple  and  inexpensive  calibration  configuration . 40 

45.  Shows  rise  time  calibration  of  two  2.082  m  long  hoses . 41 

46.  A  simple  inexpensive  diagnostic  configuration . 42 

47.  Two  short  sections  of  5/8"  permeable  hose  were  tested  using  the  end-on  test 

geometry  described  in  the  text . . 42 

48.  Simple  diagram  of  the  principle  part  of  a  two  ended  Helmholtz  oscillator . 43 

49.  A  0.32  m  section  of  5/8"  non-permeable  hose  was  set  up  as  a  Helmholtz 

oscillator  and  subjected  to  a  0.5  millibar  step  function  change  in  pressure  of 
one  manifold. . : . 44 

50.  Comparison  of  data  and  Maxhose  finite  difference  simulation  of  the 

calibration  test  configuration  of  Figure  44 . 49 


51.  Pressure  variation  in  the  hose  as  a  function  of  time  (seconds)  and  position  (X 

in  meters)  for  random  white  noise  pressure  fluctuations  in  the  atmosphere . 50 

52.  RMS  pressure  fluctuations  normalized  to  external  atmospheric  RMS 

fluctuations  at  the  end  of  the  hose  and  in  the  middle  of  the  hose  versus 
different  hose  lengths  for  a  fixed  characteristic  length  hose . 50 

53.  An  atmospheric  pressure  history  as  a  function  of  space  and  time  on  the  left 

and  the  resulting  pressure  fluctuation  signal  in  the  hose  on  the  right . 51 

54.  Pressure  fluctuations  inside  and  outside  the  middle  of  the  hose  simulation  of 

Figure  53 . 52 

List  of  Tables 

1 .  Computational  Runs . 19 


vii 


1.0  Introduction 


This  report  summarizes  work  conducted  over  a  two  year  period  to  improve  numerical 
methods.  The  report  has  three  main  sections.  Section  2  covers  work  on  synthetics  for 
layered  earth  models.  Section  3  reports  on  three-dimensional  (3D)  elastic  finite  difference 
calculations  and  the  relationship  between  one-dimensional  (ID)  (layered  earth)  and  3D 
laterally  heterogeneous  wavefields.  Section  4  reports  on  development  of  a  finite  difference 
program  to  model  the  response  of  permeable  hoses  to  atmospheric  pressure  fluctuations 
from  wind  turbulence. 

Section  2  reports  on  work  conducted  to  port  a  wavenumber  integration  synthetics 
code  to  the  Parallel  Virtual  Machine  software  and  tests  of  this  software.  Numerical 
experiments  were  performed  with  random  layering  added  to  suites  of  layered  crustal  earth 
models.  It  was  verified  that  random  layering  does  make  synthetics  look  more  like  observed 
seismograms  without  significantly  altering  the  overall  attenuation  and  spectra  of  Pg  and 
Lg  waveforms.  The  wavenumber  integration  code  was  modified  to  permit  frequency 
dependent  Q(f)  models.  Numerous  studies  have  found  that  Lg  Q(f)  is  proportional  to 
frequency  to  some  power  ranging  between  0  and  1 .  Such  studies  commonly  parameterize 
Q(f)  =  Qo  *  f1.  Numerical  experiments  were  performed  with  a  suite  of  crustal  models  to 
compare  synthetics  with  three  different  Q  models. 

In  Section  3,  we  report  on  analysis  of  the  3D  scattered  wavefield  computed  using  3D 
elastic  finite  differences  with  modal  summation  and  wavenumber  spectra  on  phase  screens. 
We  compute  scattered  waves  in  randomly  heterogeneous  3D  velocity  models.  The  total 
elastic  wavefield  (3  components  of  particle  velocity)  is  saved  on  vertical  planes  at  selected 
distances  from  the  source.  The  modal  spectra  and  wavenumber  spectra  is  estimated  on 
those  planes  and  compared  to  computations  for  a  layered  structure  in  order  to  gain  insight 
into  the  scattering  process. 

In  Section  4,  we  report  on  development  of  a  finite  difference  program,  Maxhose,  to 
compute  the  response  of  permeable  hoses  to  atmospheric  pressure  fluctuations  (noise)  and 
signals.  Permeable  hose  arrays  are  planned  as  noise  reduction  measures  for  the  deployment 
of  infrasound  stations  in  the  CTBT  International  Monitoring  System.  These  noise 
reduction  systems  are  critical  to  the  expected  performance  of  the  entire  infrasound  system. 
Maxhose  was  developed  to  provide  a  useful  tool  for  understanding  and  predicting  the 
performance  of  hoses  given  the  frequency-wavenumber  characteristics  of  turbulent  wind 
noise. 


1 


2.0  Layered  Earth  Model  Wavenumber  Integration  Synthetics 

2.1  Wavenumber  Integration  Synthetics  Using  A  Parallel  Virtual  Machine 

In  the  first  year  of  this  project  a  wavenumber  integration  program.  Prose  (Apsel  and 
Luco,  1983),  was  ported  from  a  single  processor  application  to  a  parallel  processing 
environment  using  the  Parallel  Virtual  Machine  (PVM)  software  environment.  This 
software  allows  the  user  to  harness  a  network  of  UNIX  workstations  to  perform 
calculations  in  parallel.  A  PVM  based  parallel  wavenumber  integration  program  was 
developed  that  computes  Green’s  functions  in  CSS  3.0  format  readable  with  SAC.  Near 
linear  speed-up  over  a  single  processor  has  been  realized  with  a  heterogeneous  network  of 
SUN  OS4/Solaris,  HP,  SGI,  and  DEC  workstations  tested  over  WAN  and  LAN.  The 
software  requires  PVM  3.3.11,  a  FORTRAN  F77  compiler,  and  a  C  compiler.  The  Gnu- 
make  program  is  recommended  on  SGI,  HP,  and  DEC  workstations.  Interested 
users  should  contact  Keith  L.  McLaughlin  ( scatter@maxwell.com ).  Current 
information  on  the  availability  of  the  software  can  be  found  on  the  Internet  at 
http://www.maxwell.com/products/geop.  The  interested  reader  is  referred  to  McLaughlin 
and  Shkoller  (1996)  for  detailed  descriptions  of  the  load  balancing  algorithms  employed. 

2.1.1  What  is  PVM? 

All  parallel  computer  algorithms  are  composed  of  modules  which  are  executed  on 
multiple  processors,  messages  that  must  be  sent  between  these  modules,  and  strategies  for 
coordinating  the  modules  and  processors  to  work  in  parallel.  The  modules  are  typically 
distributed  over  a  set  of  processors  that  are  connected  in  some  sort  of  communications 
network.  Much  of  the  programming  work  to  develop  a  parallel  algorithm  is  focused  on 
passing  messages  between  these  modules  and  synchronizing  their  work.  Commercially 
available  massively  parallel  processors  (MPP)  provide  custom  compilers  and  libraries  to 
facilitate  this  kind  of  programming.  However,  most  research  institutions  already  have  the 
makings  of  a  parallel  machine;  they  typically  have  several  workstations  connected  in  a 
local  area  network.  The  aggregate  computing  power  of  their  workstations  often  exceeds 
the  CPU  power  that  might  be  available  at  a  supercomputer  center.  Furthermore,  these 
machines  are  often  idle  much  of  the  day  or  night. 

Parallel  Virtual  Machine  (PVM)  is  a  public  domain  software  system  for  turning  a 
network  of  computers  into  a  virtual  parallel  computer  (Geist  et  al.  1994).  The  software 
supports  heterogeneous  networks  including  nearly  all  UNIX  workstations  and  many 
parallel  and  single  processor  supercomputers.  The  software  is  based  on  widely  used 
TCP/IP  message  passing  protocols  and  therefore  functions  over  a  wide  variety  of  local 
area  networks  (LAN)  and  wide  area  networks  (WAN).  The  user  interface  libraries  are 
both  FORTRAN  (F77  and  F90)  and  C  callable  from  a  user’s  program.  PVM  frees  the  user 
from  the  arcane  details  of  TCP/IP  message  passing  between  programs  (processes)  by 
providing  a  high  level  C  or  FORTRAN  user  interface  and  providing  buffering  and  routing 
daemons  for  to  the  TCP/IP  packets  on  the  user’s  computer.  This  is  done  by  running 
daemon  process  on  each  processor  that  serve  as  a  relay  posts  for  all  messages  passed 
between  the  user’s  programs  running  on  each  processor.  PVM  may  be  obtained  by 
anonymous  ftp  from  Oak  Ridge  National  Lab  and  University  of  Tennessee  at 
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http://www.epm.omi.gov/pviti/pvin  home.html.  Installation  of  PVM  3.3. 11. requires  an  ANSI  C 
compiler  and  can  be  accomplished  in  less  than  an  hour. 

2.1.2  Wavenumber  Integration  Synthesis  -“The  Perfect  Parallel  Algorithm” 

Wavenumber  integration  calculations  (Apsel  and  Luco,  1983)  are  computationally 
bound  by  the  time  it  takes  to  compute  the  complex  response  between  a  source  and  a 
receiver  for  a  given  frequency.  Wavenumber  integration  is  also  the  “perfect  parallel 
algorithm”  because  computation  of  one  frequency  1)  is  independent  of  all  other 
frequencies,  2)  requires  only  a  few  input  values,  3)  results  in  only  a  few  output  values,  and 
4)  there  are  many  frequencies  to  compute.  These  programs  are  well  suited  to  a  master  - 
slave  architecture  (Figure  1).  The  master  program  performs  all  file  input/output  and 
organizes  the  work  of  the  many  identical  slave  programs  running  on  multiple  processors  in 
the  network.  We  treat  each  slave  program  like  a  function  call  for  each  complex  Green’s 
function  response  at  a  fixed  frequency  and  fixed  source  and  fixed  receiver.  Figure  2 
illustrates  the  concept  behind  the  use  of  the  master  -  slave  modules  on  a  virtual  parallel 
machine  using  PVM.  The  master  program  sends  messages  to  the  slave  programs  telling 
them  which  frequencies  to  compute  and  then  waits  for  the  results  to  return  from  the  slaves 
in  the  form  of  messages.  When  the  slave  programs  are  not  computing  the  response  at  a 
specific  frequency,  they  are  waiting  for  instructions  from  the  master  program.  The  slave 
programs  return  the  numerical  results  as  well  as  timing  information  that  is  useful  in  load 
balancing  the  calculation. 
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Parallel  Virtual  Machine 


Figure  1.  The  master  program  sends  and  receives  messages  through  the  PVM 
daemon  (PVMD)  which  routes  the  messages  to  other  daemons  on  the  way 
to  and  from  slave  modules.  There  may  be  multiple  slaves  on  a  single  host 
and  there  is  often  both  a  slave  and  master  program  running  on  the  same 
host. 

Key  to  success  is  balancing  the  computational  load  between  machines  (hosts  or 
processors)  in  the  network.  A  heterogeneous  network  of  workstations  may  contain  a 
varied  mix  of  slow  and  fast  machines  Furthermore,  as  a  computation  proceeds,  the  load 
on  each  machine  will  change  with  time  as  other  users  start-up  and  terminate  other 
processes.  Also,  some  machines  may  be  on  a  local  network  and  communications  may  be 
almost  instantaneous  while  others  may  be  further  away  on  a  wide  area  network  and 
messages  may  require  a  greater  time  for  delivery.  Therefore,  it  is  necessary  to  keep  a 
concise  database  within  the  master  program  of  the  relative  performance  of  each  processor. 
We  use  a  very  simple  algorithm  for  load  balancing;  a  list  of  processors  is  ranked  by  speed 
and  a  list  of  tasks  (frequencies)  is  maintained.  Tasks  are  checked  off  the  list  as  they  are 
completed  and  the  next  task  is  always  sent  to  the  fastest  available  host.  Near  the  end  of  the 
computation,  if  we  have  sent  each  task  out  and  we  have  available  hosts,  we  do  not  wait 
for  slow  hosts  to  complete  their  tasks  but  rather  re-send  those  tasks  out  to  the  fastest 
available  host.  The  goal  is  to  keep  all  slave  modules  working.  The  resulting  algorithm  is 
therefore  robust  with  respect  to  changes  that  do  happen  in  the  network(s)  and  on  the 
various  parts  of  the  parallel  virtual  machine.  The  reader  is  referred  to  McLaughlin  and 
Shkoller  (1996)  for  details. 
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The  master  and  slave  programs,  prosem  and  proses,  are  based  on  the  wavenumber 
integration  algorithm  of  Apsel  and  Luco  (1983).  File  formats  and  program  options  are 
discussed  in  the  UNIX  style  manual  page  included  with  the  software  distribution.  The 
master  program,  prosem ,  starts  up  the  parallel  virtual  machine,  spawns  the  slave 
programs,  proses,  on  each  host  and  saves  the  frequency  domain  output  in  a  file.  A  second 
program,  gseis,  is  then  run  that  computes  the  time  domain  Green’s  functions  and  stores 
them  in  CSS  3.0  format  with  a  wfdisc  file  (see  Anderson  et  al.  1990a  and  1990b  for 
definitions  of  CSS  3. 0  format).  This  seismogram  format  can  be  read  and  manipulated  using 
a  number  of  programs  including  the  popular  Seismic  Analysis  Code  (SAC,  1995). 

2.2  Random  Layering 

It  has  been  suggested  that  introduction  of  random  velocity  variations  into  a  layered 
crustal  model  produce  more  realistic  looking  regional  phases  (Harvey,  1992).  Real 
regional  seismograms  tend  to  exhibit  extended  Pg  and  Lg  wavetrains  while  layered  Earth 
synthetics  often  have  very  simple  regional  waveforms  consisting  of  several  strong  isolated 
spikes.  One  way  to  reduce  these  isolated  spikes  and  generate  more  extended  waveforms  is 
to  introduce  interfaces  that  produce  many  small  internal  reflections  between  the  up-going 
and  down-going  waves  trapped  in  the  crust.  We  have  experimented  with  introducing  1  km 
thick  layers  and  random  velocities  distributed  about  the  mean  background  model.  The 
modified  Mooney  et  al.  (1997)  layered  Earth  models,  and  details  of  the  synthetic 
computation  procedures  are  described  in  Bennett  et  al.  (1997),  and  the  frequency  and 
depth  dependent  attenuation  model,  Q(f,z)  =  Q0(J3,z)f1,  is  described  in  Section  2.3  of  this 
report.  Figures  2  and  3  compare  Green's  functions  model  D2  with  and  without  random 
layers  of  5%  variation.  The  Pn  and  Sn  waveforms  in  the  two  sets  of  Green's  functions  are 
larger  and  more  impulsive  with  the  random  layering  than  without.  Also,  some  of  the 
isolated  spikes  in  the  Pg  and  Lg  wavetrains  are  reduced. 

Figure  4  shows  more  such  random  perturbations  to  the  DA,  GA,  H2,  and  N4  models 
of  Mooney  et  al.  (1997).  Variations  of  0%,  2.5%,  5%,  7.5%,  and  10.0%  RMS  random 
velocity  variation  have  been  introduced  into  both  P  and  S  wave  velocities  throughout  each 
model.  Green's  functions  for  these  models  are  compared  in  Figure  5  through  Figure  8.  The 
four  background  models  were  chosen  for  their  diversity;  the  models  represent  a  thick  set 
of  platform  deposits  (DA),  a  high-velocity  crust  with  very  thin  sediments  (GA),  a 
moderate  velocity  crust  with  no  sediments  (H2),  and  a  thin  crust  with  sediments  (N4). 
Smoothed  Lg/Pg  spectral  ratios  computed  from  these  synthetic  seismograms  are  shown  in 
Figure  11. 

It  is  ironic  that  while  our  principal  goal  was  to  generate  more  realistic  Pg  and  Lg 
waveforms,  introduction  of  random  layering  has  made  the  Pn  and  Sn  waveforms  larger, 
more  impulsive,  and  more  realistic.  These  waves  travel  in  the  upper  mantle  waveguide  and 
are  sensitive  to  the  velocity  gradients  just  below  the  Moho.  The  background  mantle  PEM 
model  has  no  gradient  except  for  the  “earth  flattening”  approximation.  We  presume  that  in 
the  absence  of  the  many  interfaces,  the  waves  consist  of  only  a  few  rays  or  modes  that 
propagate  in  the  upper  mantle  waveguide.  Addition  of  the  random  interfaces  appears  to 
proliferate  the  number  of  rays/modes  that  contribute  to  the  Pn  and  Sn  waveforms. 
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Uxy  Greens  Function,  h  =  15  km,  Delta  =  400  km 


Seconds 


Figure  2.  Comparison  of  synthetic  transverse  component  Green's  function,  Gyxy,  at 
400  km,  source  depth  h=15km,  with  (top)  and  without  (below)  5%  random 
velocity  variations  in  the  layered  model  D2.  Note  that  the  more  impulsive  Sn 
arrival  from  the  random  layered  structure  and  that  the  reduced  isolated 
"spikes"  in  the  Lg  signal. 


Explosion  Greens  Function,  h  =  15  km,  Delia  =  400  km 


100  120 
Seconds 


Figure  3.  Comparison  of  vertical  component  Green's  function,  G*,  at  400  km,  source 
depth  h=lkm,  with  (top)  and  without  (below)  5%  random  velocity 
variations  in  the  layered  model  D2.  Note  the  more  impulsive  Pn  waveform 
and  the  reduced  "spikes"  in  the  Pg  from  the  randomized  structure. 
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Model  DA  wilh  0. 2.5,  5.0  7.5  end  10%  RMS  ModolGA  wilh  0. 2.5. 5.0  7.5  and  10%  RMS 


Figure  4.  P-wave  velocities  versus  depth  for  layered  models  DA,  GA,  H2,  and  N4 
with  0%,  5%,  7.5%,  and  10%  velocity  variations. 
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Comparison  of  Green's  functions  for  model  DA  with  0%,  2.5%,  5.0%, 
7.5%,  and  10%  velocity  variation  (top  to  bottom). 
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Comparison  of  Green's  functions  for  model  GA  with  0%,  2.5%,  5.0%, 
7.5%,  and  10%  velocity  variation  (top  to  bottom). 


Comparison  of  Green's  functions  for  model  H2  with  0%,  2.5%,  5.0%,  7.5%, 
and  10%  velocity  variation  (top  to  bottom). 


Comparison  of  Green's  functions  for  model  N4  with  0%,  2.5%,  5.0%,  7.5%, 
and  10%  velocity  variation  (top  to  bottom). 
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Figure  11.  Vertical  component  Lg/Pg  spectral  ratios  for  models  DA,  GA,  H2,  and  N4 
for  an  explosion  source  (h=lkm)  and  0, 2.5,  5,  7.5,  and  10%  RMS  random 
velocity  layers.  Note  that  the  general  character  of  the  Lg/Pg  ratio  while  very 
different  for  each  background  model,  is  not  altered  much  by  the  introduction 
of  random  velocity  layers. 

2.3  Frequency  Dependent  Attenuation,  Q(f)  =  Qo*r 

We  have  examined  three  simple  crustal  attenuation  models,  Ql,  Q2,  and  Q3  in  detail. 

Ql)  Qn(f,z)  =  Qo/f1, r|  =  0,  Qon  =  (3(z)/10  for  z  >  0  m, 

Q2)  Qn(f,z)  =  Qo/f1,  ti  =  0,  Qon  =  P(z)/5  for  z  >  3  km,  low  Q  surface  layers 
Q3)  CM£z)  =  Qo/f1,  ri  =  0.5,  Q0r  =  P(z)/5  for  z  >  3  km,  low  Q  surface  layers. 

The  choice  of  model  Q3  was  based  on  three  requirements. 

•  Low  Q  surface  layers  are  required  to  attenuate  unwanted  short-period  low-group 
velocity  higher  modes  and  Rg  that  propagate  in  the  low-velocity  near  surface 
layers.  We  chose  Q0m.  =  25  for  0  <  z  <  1  km,  Q0h  =  75  for  1  km  <  z  <  2  km  and  Q0(l 
=  150  for  2  km  <  z  <  3  km. 

•  Average  shear  and  compressional  Q  values  in  the  crust  must  generally  produce 
Lg/Pg  ratios  greater  than  unity  near  1  Hz  for  earthquake  mechanisms.  We  found 
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that  Qon  =  P(z)/5  for  z  >  3  km  would  yield  reasonable  Lg/Pg  ratios  near  1  Hz, 
while  the  model  Q0fl  =  0(z)/lO  attenuated  Lg  too  much. 

•  Short-period  earthquake  Lg/Pg  spectral  ratios  as  a  function  of  frequency  should 
generally  remain  above  or  near  unity  from  1  to  5  Hz.  We  chose  Q  proportional  to 
f0'5  in  order  to  keep  Lg/Pg  ratios  from  declining  too  steeply  as  a  function  of 
increasing  frequency. 

Figures  12  and  13  show  compilations  of  Lg/Pg  spectral  ratios  from  Bennett  et  al. 
(1997).  Earthquake  Lg/Pg  ratios  generally  slowly  decrease  with  increasing  frequency. 

Synthetics  from  a  variety  of  crustal  models  suggested  that  in  order  to  produce 
earthquake-like  Lg/Pg  ratios  greater  than  unity  near  1  Hz  the  higher  Qo’s  of  models  Q2 
and  Q3  are  preferred  to  those  of  model  Ql.  Also,  to  keep  Lg/Pg  ratios  relatively  flat  as  a 
function  of  frequency,  an  increase  in  Q  as  a  function  of  frequency  (n  >  0)  is  needed.  The 
Lg/Pg  ratios  above  1  Hz  for  models  Ql  and  Q2  were  too  small  and  do  not  agree  with 
observations.  Lg  Q  increasing  with  increasing  frequency  is  commonly  observed  (Nuttli, 
1981;  Goncz  et  al,  1986,  Gupta  and  McLaughlin  1987;  Campillo  et  al.,  1985;  Mitchell 
1981)  with  r\  between  0  and  1.  It  should  be  noted  that  numerous  researchers  have  found  a 
negative  correlation  between  Q0  and  r\;  the  higher  Q  is  at  1  Hz,  the  slower  it  increases 
with  increasing  frequency.  Therefore,  it  may  be  possible  to  reproduce  many  of  the 
observed  results  by  assuming  higher  Q0's  with  somewhat  smaller  values  of  r|  <  0.5.  Also, 
in  tectonic  regions  with  lower  1  Hz  Q  values,  the  value  of  t|  may  be  higher. 

Figures  15-17  show  explosion  and  earthquake  Green's  functions  computed  for  the 
three  Q  models.  The  velocity  model  D2  is  well  suited  to  show  off  the  differences  between 
Ql,  Q2,  and  Q3.  The  earthquake  Lg  amplitudes  are  too  small  in  models  Ql  and  Q2  for 
both  the  broadband  and  the  higher  frequency  seismogram  (0.5  Hz  highpass  Figure  14  and 
3.0  Hz  highpass  Figure  16).  A  shallow  explosion  excites  very  slow  unrealistic  waves  in  the 
near  surface  layers  in  model  Ql  (Figure  15). 
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Compilation  of  average  explosion  Lg/Pg  spectral  ratios  from  Bennett  et  al. 
(1997)  for  explosions  in  Asia  and  Western  US.  Note  that  the  average  Lg/Pg 
is  generally  greater  than  unity  for  frequencies  at  and  below  1  Hz. 


Average  Earthquake  Lg/Pg  Ratios  from  Bennett  et  al.  (1997) 
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Figure  13.  Compilation  of  average  Lg/Pg  spectral  ratios  from  Bennett  et  al.  (1997)  for 
earthquakes  in  Asia,  Western  US,  and  Eastern  US.  Note  that  Lg/Pg  is 
generally  greater  than  unity  for  frequencies  below  5  Hz  with  the  exception 
of  the  ARU  data.  The  average  earthquake  curves  are  generally  above  the 
corresponding  explosion  curves  at  frequencies  above  1  Hz. 
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Figure  14.  Comparison  of  vertical  component  synthetics  (high  pass  at  0.5  Hz)  at  200 
km  for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle), 
Q3  (bottom).  A  dip-slip  double-couple  source  (Gzds)  at  15  km  depth  is 
shown.  The  Lg  amplitude  is  significantly  increased  relative  to  the  Pg 
amplitude  for  t|=0.5,  Qo=|3/5,  model  Q3. 
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Figure  15.  Comparison  of  vertical  component  synthetics  (high  pass  at  0.5  Hz)  at  200  km 

for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle),  Q3 
(bottom).  An  explosive  source  (Gzi)  at  1  km  depth  is  shown.  The  high  Q 
surface  layer  model  does  not  sufficiently  attenuate  the  late  arriving 
fundamental  surface  waves  excited  by  the  shallow  explosive  source.  Also,  it 
is  clear  that  the  Ql,  high  Q  surface  layer  model  contains  too  much  shallow 
propagating  energy  in  the  Lg  window  that  is  absent  in  the  low  Q  surface 
layer  models,  Q2  and  Q3. 
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Figure  16.  Comparison  of  vertical  component  synthetics  (high  pass  at  3.0  Hz)  at  200 
km  for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle), 
Q3  (bottom).  A  dip-slip  double-couple  source  (Gzds)  at  15  km  depth  is 
shown.  The  Lg  amplitude  is  significantly  increased  relative  to  the  Pg 
amplitude  for  r|=0.5,  Qo=p/5,  model  Q3 


Figure  17.  Comparison  of  vertical  component  synthetics  (high  pass  at  3.0  Hz)  at  200 
km  for  three  Q(f,z)  models  and  velocity  model  D2,  Q1  (top),  Q2  (middle), 
Q3  (bottom).  An  explosive  source  (Gzi)  at  1  km  depth  is  shown  below.  The 
Lg  is  almost  non  existent  in  model  Q2. 
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3.0  Finite  Difference  3D  Regional  Scattering  Calculations 

3.1  Introduction 

Numerical  methods  for  3D  elastic  wave  propagation  have  made  significant  advances  in 
the  last  decade  due  in  part  to  better  software,  but  largely  due  to  exponential  growth  of 
computer  capabilities.  For  a  fixed  price,  computer  memory  and  speed  roughly  double 
every  18  months  (popularly  known  as  "Moore's  Law").  Provided  such  technological 
advances  continue  it  will  be  another  decade  or  more  before  3D  computations  are  as 
economical  as  ID  computations  today.  3D  computations  require  8  times  more  memory 
and  16  times  more  CPU  time  than  comparable  ID  computations.  For  example, 
computation  by  wavenumber  integration  of  a  record  section  of  0.5  Hz  regional  Green's 
functions  from  0  to  200  km  at  0.5  km  intervals  requires  about  l/20th  the  time  as  a  section 
of  3D  waveforms  for  a  single  source  (hours  versus  days).  While  the  suite  of  3D 
seismograms  is  far  richer  than  the  ID  record  section,  the  bandwidth  and  range  of  the  3D 
computation  is  far  more  limited  than  the  ID  computation.  With  this  in  mind,  it  is  beneficial 
to  leverage  ID  numerical  methods  as  much  as  possible  in  understanding  3D  results  and 
devising  hybrid  methods  for  seismogram  calculation. 

Seismologists  have  developed  extensive  insight  into  the  nature  of  wave  propagation  in 
layered  earth  models  using  ID  methods.  Layered  earth  model  synthetics  correctly  predict 
many  of  the  gross  features  of  regional  seismograms.  Distinct  regional  phases  such  as  Pn, 
Pg,  Sn,  Lg,  and  Rg  propagate  as  loosely  coherent  wave  packets  with  well  defined  group 
velocities  for  large  distances  in  crustal  and  upper  mantle  waveguides.  3D  methods  aim  to 
predict  how  these  phases  are  affected  by  scattering.  For  some  cases  we  can  view  the  3D 
scattering  as  a  perturbation  upon  the  waveguides.  From  the  point  of  view  of  modal 
summation  methods,  we  can  examine  how  the  "modes"  defined  by  the  average  ID  layered 
structure  exchange  energy  with  increasing  range.  Alternatively,  we  can  view  the  outgoing 
wavefield  from  the  point  of  view  of  wavenumber  integration  or  phase  screen  methods.  We 
can  view  the  scattering  as  converting  energy  from  one  wavenumber  into  energy  at  a 
different  wavenumber. 

In  this  section,  we  report  on  work  we  have  done  to  analyze  the  scattered  wavefield 
computed  using  3D  elastic  finite  differences  with  modal  summation  and  wavenumber 
spectra  on  phase  screens.  We  compute  scattered  waves  in  3D  velocity  models  constructed 
by  adding  random  perturbations  to  a  mean  crustal  model.  The  total  elastic  wavefield  (3 
components  of  particle  velocity)  is  saved  on  vertical  planes  at  selected  distances  from  the 
source.  We  then  estimate  the  modal  spectra  and  the  wavenumber  spectra  on  those  planes. 
We  compare  these  spectra  with  computations  for  a  layered  structure  in  order  to  gain 
insight  into  the  scattering  process. 

Tres3D  with  recursive  grid  refinement  was  used  to  compute  complete  3D  elastic  wave 
propagation  in  a  3D  laterally  heterogeneous  model.  Tres3D  is  a  2nd  order  explicit  finite 
difference  code  for  3D  rectangular  meshes.  It  requires  a  minimum  of  10  cells/wavelength 
and  the  grids  were  designed  to  provide  a  bandwidth  from  0  to  0.6  Hz.  Sources  are 
inserted  as  moment  tensors.  Recursive  grid  refinement  (RGR)  described  in  McLaughlin 
and  Day  (1995)  uses  nested  grids  to  decrease  memory  and  CPU  time.  Fine  zoning  is  used 


17 


in  the  crust  (250  m),  coarser  zoning  is  used  in  the  upper  most  mantle  (500  m)  and  the 
coarsest  zoning  (1000  m)  is  used  in  the  deepest  mantle.  Moderate  to  high  attenuation  is 
used  in  the  mantle  to  suppress  reflections  from  the  bottom  of  the  grid  and  the  coarser 
grids.  The  sources  were  placed  at  X  =  0,  Y  =  0  at  depths  of  Z  =  -0.5  km  and  Z  =  -12.5 
km.  A  reflection  symmetry  axis  was  used  as  a  boundary  condition  on  the  Y  =  0  and  X  =  0 
planes  to  reduce  the  size  of  the  problem  by  a  factor  of  4.  The  geometry  is  diagrammed  in 
Figure  3.1-1. 

The  computation  was  performed  over  a  200  by  100  by  100  km  volume  and  run  to 
durations  of  about  100  seconds.  The  recursive  grid  refinement  was  used  with  21  grids  and 
a  total  of  10  million  cells.  Each  computation  required  about  3  CPU  days  on  a  DEC  2100 
5/250  workstation  to  complete  the  7.3*10n  cell-cycles. 

A  list  of  computations  is  given  in  Table  1.  Two  source  types  were  modeled:  an 
explosion  source  with  Mxx  =  Myy  =  Mzz  centered  at  a  depth  of  0.5  km  and  a  "double¬ 
couple"  source  with  Mxx  =  -Mzz,  Myy  =  0  at  a  depth  of  12.5  km.  Lateral  heterogeneity 
was  limited  to  the  crust  and  upper  most  mantle.  The  mantle  was  a  homogeneous  half¬ 
space  and  a  simple  linear  gradient  was  chosen  for  the  background  crustal  structure.  The 
background  velocity  model  is  shown  in  Figure  18.  In  order  to  introduce  lateral 
heterogeneity  we  used  the  following  procedure: 

1)  Generate  a  3D  random  array,  r(x,y,z),  uniformly  distributed  between  -1  and  1, 
sampled  at  x  =  0,  ....w,  y  =  0,  ...,  ymax,  z  =  z™,  ...,0  to  fill  the  volume  xmax  = 
128  km,  y^  =  32  km,  z™,  =  -32  km,  with  sampling  intervals  dx  =  1.0  km  ,  dy  = 
1.0  km,  and  dz  =  0.5  km.. 

2)  Smooth  the  random  numbers  in  the  X  and  Y  directions  with  a  3-point  smoothing 
operator  (1,1,1).  Re-scale  the  array  to  an  RMS  value  of  1. 

3)  For  a  given  (x,y,z)  location  in  the  crust  (z  >  Zmin)  set  the  S-velocity  to  p  x(x,y,z)  = 
P_1o(z)  (  1  +  RMS  *  r(x,y,z)  ).  For  x  >  x^,  or  y  >  y^,  or  z  >  z^,  the  velocity  is 
set  to  the  background  model.  Run  #1  was  formulated  with  fluctuations  specified  in 
velocity  instead  of  slowness. 

This  procedure  introduces  an-isotropic  heterogeneity  with  characteristic  lengths  of  about  1 
km  in  the  Z  and  X  directions  and  2  km  in  the  Y  direction.  Figures  20  and  21  show  two 
ways  to  visualize  this  lateral  heterogeneity. 

Snap  shots  of  vertical  particle  velocity  from  Run#l  (5%  RMS  velocity  variation, 
explosion  source)  along  the  Y  =  0  plane  are  shown  in  Figures  22  through  27.  Complexity 
of  the  wave  propagation  is  immediately  evident.  Development  of  the  refracted  P  and  Pn 
can  be  seen  at  T  =  15  seconds  as  well  as  development  of  the  moho  reflection,  PmP.  Lg  can 
be  best  seen  at  later  times  T  =  20,  25,  and  30  seconds  as  an  interference  of  up-going  and 
down-going  waves  in  the  crust.  The  short-period  fundamental  Rayleigh  wave,  Rg,  can  be 
seen  running  along  the  surface  at  about  2.5  km/s.  The  disorganized  wavefield  behind  the 
Rg  wavefront  is  simply  labeled  "coda".  Waveforms  from  Run#l  are  shown  in  Figure  28  to 
illustrate  that  Rg  is  the  most  significant  arrival  at  these  distances  on  the  surface,  and  it  is 
systematically  delayed  and  attenuated  by  the  random  fluctuations. 
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Table  1.  Computational  Runs* 


Run# 

Source  Type 

RMS  Variation  : 

Source  Dei 

Run  1 

Mxx  =  Myy  =  Mzz 

5%  velocity 

0.25 

Run  2 

Mxx  =  Myy  =  Mzz 

5%  slowness 

0.25 

Run  3 

Mxx  =  Myy  =  Mzz 

5%  slowness 

0.25 

Run  4 

Mxx  =  Myy  =  Mzz 

0 

0.25 

Run  5 

Mxx  =  Myy  =  Mzz 

5%  slowness 

0.25 

Run  6 

Mxx  =  Myy  =  Mzz 

7.5%  slowness 

0.25 

Run  7 

Mxx  =  -Mzz,  Myy  =  0 

7.5%  slowness 

12.5 

Run  8 

Mxx  =  -Mzz,  Myy  =  0 

0 

12.5 

200  km 


■■■•>•  .s  \ 


•100  km 


''?S'<fx=O.Skm 


S-wave  velocities 
‘  0.25.  The  Moho 


Visualization  of  the  3D  random  variation.  Isosurfaces  with  absolute 
variation  greater  ±lc  are  shaded  in  the  top  Image.  Regions  within  ±ler  are 
transparent.  The  lower  image  shows  the  same  with  ±3a.  These  images 
help  visualize  the  cores  of  the  scatterers.  Note  the  elongation  of  the  scatter 
cores  in  the  Y  direction. 


Vertical  velocity,  Y=©  plane,  T=5  sec,  5%  RMS 
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Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  5 
seconds. 


Vertical  velocity,  Y=0  plane,  T=10  sec,  5%  RMS 

0 

-10000 

£ 

^  -20000 
-30000 

0  25000  50000  75000  100000 

X(m) 

Figure  23.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  10 

seconds. 


Vertical  velocity,  Y=0  plane,  T=15  sec,  5%  RMS 


Figure  24. 
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Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  15 
seconds. 


Vertical  velocity,  Y=0  plane,  T=20  sec,  5%  RMS 
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Figure  25.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  20 
seconds. 


Vertical  velocity,  Y=0  plane,  T=25  sec,  5%  RMS 
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Figure  26.  Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  =  25 
seconds. 

Vertical  velocity,  Y=G  plane,  T=30  sec,  5%  RMS 
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Figure  27. 
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Run  #1  snap  shots  of  vertical  particle  velocity  on  the  Y  =  0  plane,  at  T  = 
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Figure  28.  Comparison  of  finite  difference  and  wavenumber  integration  synthetics  for 
Run  #1  (5%  RMS,  explosion  source)  at  100  km  for  Iowpass  filters  0.2,  0.4, 
0.6,  and  0.8  Hz  (top  to  bottom).  Note  that  the  Rg  is  unaffected  in  the  lowest 
bandwidth  and  is  increasingly  attenuated  and  delayed  with  increasing 
frequency.  Frequencies  above  0.6  Hz  are  affected  by  grid  dispersion. 

3.2  Modal  Spectra 

Modal  summation  is  a  common  approach  to  either  analyze  or  compute  wave 
propagation  in  a  layered  structure  (Aki  and  Richards,  1980).  The  methods  for  seismogram 
synthesis  have  a  long  and  fruitful  history.  In  particular,  Lg  can  be  analyzed  as  groups  of 
the  P-SV  and/or  SH  modes.  Several  researchers  have  described  weak  scattering  as  mode¬ 
mode  conversion,  where  energy  is  transferred  from  one  mode  to  another.  We  wish  to  take 
advantage  of  the  machinery  developed  for  modal  summation  and  estimate  a  "modal 
spectra"  or  a  "modal  decomposition".  We  write  the  vertical  seismogram  at  frequency,  f, 
depth,  z,  and  distance.  A,  as  a  mode  sum,  Sz(f,z,A)  =  Sj=j,N  Czj(f,A)  Ezj(z,f),  where  Ezj(z,f) 
is  the  j'th  mode  vertical  eigenfunction  at  frequency,  f  (see  Figures  29  and  30).  Likewise 
for  the  radial  seismogram,  Sa(f,z,A)  =  2j=j, ...  ,n  CRj(f,A)  Eaj(z,f).  We  can  then  examine  the 
modal  spectra,  Cj(f,A),  as  a  function  of  frequency  and  distance.  In  order  to  estimate  this 
modal  spectra,  we  save  the  vertical  and  radial  velocities  at  some  distance  for  each  depth 
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in  the  finite  difference  grid.  Since  the  eigenfunctions  are  orthogonal,  we  should  only  need 
to  compute  an  inner  product  between  S(f,z,A)  and  Ej(z.f)  to  compute  the  coefficients, 
Cj(f,A).  However,  the  finite  difference  seismograms  are  only  sampled  at  specific  depths 
from  the  surface  to  a  maximum  depth  and  they  contain  wave  types  that  are  not  present  in 
the  pure  P-SV  modes,  such  as  Pn  and  Pg  as  well  as  scattered  body  waves.  In  short,  the 
modes  are  not  exactly  orthogonal  under  the  numerical  quadrature  and  leakage  can.  occur 
of  the  non  modal  waves  into  the  modes.  It  helps  to  use  both  the  vertical  and  radial  motion 
to  constrain  the  modal  estimates.  We  use  a  stripping  procedure  described  below. 

1)  start  with  mode]  =  1  and  set 

S’z(f,z,A)  =  Sz(f,z,A) 

S’R(f>z,A)  =  SR(f,z,A) 

2)  estimate  the  normalized  inner  product  over  depth  from  z\  to  znz 

EzjoCO  =  Ek=i,nz  Ezj(Zk;f)“ 

ERjo(f)  =  Ek=RRZ  ERj(Zk,f}~ 

Cj(f,A)=(l/2)Sk=3.nz(S’z(f,zk,A)Ezj(zk,f),/Ezjo(f)+S’R(fJZk,A)ERj(zk,f)/ERj0(f))dzk 

3)  strip  the  j'th  component  from  the  observed  seismogram 

S'z(f,z,A)  =  S'z(f,z,A)  -  Cj(f,A)  Ezk(z,f) 

S'R(f,z,A)  =  S’R(f,z,A)  -  Cj(fiA)  ERj(z,f) 

4)  if  j  <  N  then  increment  the  mode  number]  to  j  +  1  and  go  to  step  2. 

We  find  this  stripping  procedure  is  relatively  insensitive  to  whether  we  start  with 
mode  j  =  1  and  increment  to  mode  j  =  N,  or  start  with  mode  j  =  N  and  decrement  to  mode 
j  =  1.  Amplitude  of  the  residual  seismogram  is  small  in  comparison  to  the  original 
seismogram.  Generally,  about  90%  of  the  seismogram  energy  is  accounted  for  by  the 
modal  spectra.  Obviously,  this  simple  quadrature  rule  could  be  refined,  but  we  believe 
that  the  modal  spectra  estimated  in  this  way  provide  insight  into  the  scattering  processes. 
Figure  3 1  shows  modal  spectra,  Cj(f,A),  for  Run#2  (5%  RMS)  and  Run#4  (0%  RMS)  at  a 
range  of  120  km.  Note  the  modal  cut-off  where  the  number  of  modes  increase  with 
increasing  frequency.  Mode  1  corresponds  to  the  fundamental  Rayleigh  wave.  Modes  2 
through  8  are  the  higher  modes  that  exist  in  the  test  structure  for  frequencies  up  to  0.6  Hz. 
The  modes  are  normalized  to  unit  vertical  amplitude  at  the  free  surface.  Therefore,  the 
spectral  amplitudes  can  be  interpreted  in  terms  of  the  amplitudes  of  the  modes  of  an 
observed  surface  vertical  seismogram.  The  modes,  however,  may  not  necessarily  arrive  at 
their  customary  group  velocities  nor  as  isolated  modes  in  the  time  domain.  Clearly,  the 
fundamental  is  the  largest  amplitude  mode  from  0.1  to  about  0.4  Hz  for  both  the  laterally 
homogeneous  and  heterogeneous  models.  The  fundamental  is  the  largest  mode  for  the 
laterally  homogeneous  model  (Run#4)  at  all  frequencies,  but  has  been  reduced  to  an 
amplitude  comparable  to  the  higher  modes  at  frequencies  greater  than  about  0.4  Hz  in  the 
laterally  heterogeneous  model  (Run#2).  We  interpret  this  as  scattering  from  the 
fundamental  into  higher  modes  by  lateral  heterogeneity.  A  simple  way  to  illustrate  this  is 
shown  in  Figure  32  where  we  form  a  ratio  between  the  modal  spectra  from  Run#2  and 
Run#4.  This  modal  spectral  ratio  demonstrates  the  enhancement  of  energy  in  the  higher 
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modes  at  nearly  all  frequencies  at  the  expense  of  the  fundamental.  Likewise,  we  show  a 
ratio  of  modal  spectra  at  60  km  to  spectra  at  120  km  in  Figure  33. 
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Figure  29.  Plots'  of  P-SV  modes  at  0.5  Hz  (right)  and  0.75  Hz  (left)  for  the  test 
structure  shown  in  Figure  19.  Note  that  the  fundamental  mode  has  very  little 
motion  below  4  km  (at  0.75  Hz)  and  below  5  km  (at  0.5  Hz).  The  higher 
modes  sample  the  entire  crust. 


g 

s; 


8 


(J: 


Mode  N 


Mode 


Cut 


umber,  J 


KM 


Figure  30.  Diagram  of  Modal  Spectra.  C(J,F).  C(T,F)  =  0  to  the  right  of  the  modal  cut¬ 
offline. 
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Modal  Spectra,  120  km  Vertical  Component,  0%  RMS  Modal  Spectra,  120  to  Vertical  Component,  5%  RMS 
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Modal  Spectra,  60  to  Vertical  Component,  0%  RMS  Modal  Spectra,  60  to  Vertical  Component,  5%  RMS 
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Figure  31.  Modal  spectra  at  120  km  (top)  and  60  km  (bottom)  for  models  with  0% 
RMS  (left)  and  5%  RMS  (right)  lateral  heterogeneity. 
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Ratio  Modal  Spectra,  120  km,  5%  RMS/0%  RMS 
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Figure  32.  Modal  Spectral  Ratio  at  120  km  5%  RMS  model  /  0%  RMS  model.  Yellow 
and  Red  shades  are  modes  enhanced  w.r.t.  the  background  model.  Blue  and 
purple  shades  are  modes  deficient  w.r.t.  the  background  model. 


Ratio  Modal  Spectra,  120  km  to  60  km 
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Figure  33.  Ratio  of  120  km  to  60  km.  This  spectral  ratio  highlights  modes  that  have 
lost  or  gained  energy  between  60  and  120  km. 
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3.3  Wavenumber  Spectra 

While  the  modal  spectra  are  useful  descriptions  of  the  crustal  wavefield,  they  do  not 
provide  a  complete  description.  In  order  to  obtain  this  more  complete  description  of  the 
wavefield,  we  have  taken  a  cue  from  other  numerical  methods  for  computing  wavefields 
such  as  wavenumber  integration  or  phase  screen  methods.  We  save  the  3  components  of 
motion  on  a  32  km  by  32  km  vertical  plane  perpendicular  to  the  nominal  direction  of 
propagation  (X  =  constant)  and  compute  a  3D  FFT  of  the  seismograms  to  transform  them 
from  the  (y,  z,  i)  domain  to  the  wavenumber-frequency  (Fy  ,FZ,  f)  domain.  Fz  =  f  sz,  and  Fy 
=  f  Sv  are  the  wavenumbers  for  propagation  in  the  vertical  (Z)  and  transverse  (Y) 
directions  respectively,  sz  and  Sy  are  the  slowness  components  in  the  vertical  and 
transverse  directions  respectively. 


S(Fy,Fz,f)  = 


The  diagram,  in  Figure  34  can  be  of  help  in  understanding  these  wavenumber- 
wavenumber  spectra.  Note  that  because  the  finite  difference  calculations  are  performed 
with  a  symmetry  axis  at  Y  -  0,  the  phase  plane  spectra  are  also  symmetric  about  the  Z 
axis,  S(Fv,  Fr)  =  S(-Fy,  Fz).  Waves  traveling  straight  at  the  phase  plane  are  plotted  at  the 
origin  of  the  diagram.  Upgoing  waves  arriving  perpendicular  to  the  plane  plot  along  the 
positive  Fz  axis,  while  downgoing  waves  arriving  perpendicular  to  the  phase  plane  plot 
along  the  negative  Fz  axis.  Off  azimuth  waves  (side-swipe)  plot  away  from  the  Fz  axis. 
Since  there  is  a  minimum  velocity  in  the  crustal  waveguide,  all  propagating  waves  must 

plot  inside  a  circle  given  by  (Fz2  +  Fy2)in  =  f  /  {U,.  Waves  with  Fz  =  f  s*  <  frP-oniie  are 

trapped  in  the  crust  where  Pmontie  ls  the  shear  wave  velocity  of  the  mantle  below  the 
moho.  Because  the  phase  plane  has  finite  dimensions,  there  are  side  lobes.  Also,  all 
outgoing  waves  are  not  planar  and  they  have  curvature  as  they  impinge  upon  the  phase 
plane;  even  in  the  absence  of  lateral  heterogeneity’,  the  spectra  are  not  confined  strictly  to 
the  Fz  axis.  Therefore,  for  each  wavenumber  spectra  for  a  laterally  heterogeneous  model 
we  show  the  corresponding  wavenumber  spectra  for  the  strictly  layered  structure  for 
comparison.  And,  in  order  to  gain  insight  into  the  side  lobes  and  impulse  response  to  a 
curved  wavefront  at  a  selected  distance,  we  plot  the  impulse  response  of  some  cylindrical 
waves  incident  upon  the  32  km  by  32  km  phase  plane  at  120  km  from  a  point  source  in 
Figures  35  and  36. 
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Figure  34.  Diagram  of  wavenumber-wavenumber  spectra  at  frequency,  f.  All 
wavenumber  spectra  in  this  paper  are  symmetric  about  the  vertical 
wavenumber  axis  because  of  a  reflection  symmetry  on  the  Y  =  0  plane  in 
the  finite  difference  calculations. 
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Figure  35.  Phase  plane  wavenumber  impulse  response  for  a  cylindrical  wave,  120  km 
from  the  source  at  0.4  (left)  and  0.6  Hz  (right). 


mm 


Bange=120  km,  FrO£Hz>  Off  Axis  Cylindrical  Wave 


N 

u. 


Fy  (1/m) 


~60. 0~5G. 0-40  i)"3D.O"20i>~1 0,0 


20*log(U) 

Figure  36.  Phase  plane  wavenumber  impulse  response  for  an  off-axis  cylindrical  wave, 
120  km  from  the  source  at  0.6  Hz. 

Wavenumber-wavenumber  spectra  at  120  km  and  0.6  Hz  from  Run  #  2  (5%  RMS 
explosive  source)  and  Run  #4  (0%  RMS,  explosive  source)  are  compared  in  Figures  37- 
39.  In  the  absence  of  scattering,  on-azimuth  up-going  and  down-going  waves  with 
vertical  slownesses  near  0.3  sec/km  (apparent  vertical  velocity  about  3  km/sec)  are 
particularly  strong  on  the  vertical  and  radial  components  of  motion  (Z  and  X). 
Presumably  these  correspond  to  the  developing  Lg  wavetrain.  Energy  with  vertical 
slownesses  less  than  0.2  sec/km  is  presumably  related  to  the  developing  Pg  wavetrain  as 
well  as  steeply  incident  P  and  S  waves  that  should  escape  into  the  mantle.  For  all  three 
components  of  motion,  we  see  that  scattering  has  stretched  out  the  spectra  in  the  Fy 
direction  indicative  of  the  off-azimuth  energy.  Scattering  has  generally  lowered  the  on- 
azimuth  energy.  Note  that  the  transverse  particle  motion  (Y  component)  is  not  identically 
zero  for  the  laterally  homogeneous  case  since  the  X  =  120  km  plane  is  only  perpendicular 
to  the  direct  waves  at  Y  =  0.  The  scattered  Y  component  wavenumber  spectra  become 
very  complicated  and  fill  a  larger  wavenumber  region. 
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Figure  38.  Wavenumber  spectra  of  the  radial  component  velocity  field  on  a  Y  =  120  km 
phase  plane  at  0.6  Hz  for  the  5%  RMS  (right)  and  0%  RMS  (left). 
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Figure  39.  Wavenumber  spectra  of  the  transverse  component  velocity  field  on  a 
Y  =  120  km  phase  plane  at  0.6  Hz  for  the  5%  RMS  (right)  and  0%  RMS 
(left). 

3.4  Conclusions 

Modal  and  wavenumber  spectra  of  phase  planes  provide  tools  to  analyze  finite 
difference  simulations  of  regional  scattering.  Modes  serve  as  a  convenient  set  of  basis 
functions  and  have  a  straight  forward  interpretation.  Modal  spectral  ratios  can  be  used  to 
estimate  scattering  attenuation  of  each  mode  and  the  general  transfer  between  modes. 
Wavenumber  spectra  are  perhaps  more  difficult  to  visualize  but  they  provide  means  to 
visualize  off-azimuth  waves  and  waves  that  are  not  necessarily  trapped  in  the  crustal 
waveguide. 

Fundamental  mode  dispersion  looks  significantly  different  for  the  randomized  models 
than  predicted  by  the  average  model.  The  scattering  both  attenuates  the  Rg  and  introduces 
significant  group  delay.  If  scattering  makes  a  significant  contribution  to  fundamental 
Rayleigh  attenuation,  seismologists  should  be  careful  interpreting  the  apparent  group 
velocities  and  the  inferred  shallow  velocities. 

Conversion  from  the  fundamental  Rg  becomes  increasingly  important  with  increasing 
frequency.  With  only  moderate  heterogeneity  (5%  RMS)  the  P-SV  modes  are  rapidly 
approaching  equilibrium  at  distances  less  than  100  km  for  frequencies  above  0.5  Hz. 
Likewise,  both  off-azimuth  P-SV  and  SH  energy  is  significant  at  60  km. 
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4.0  Permeable  Hose  Characteristics  and  Noise  Reduction 
For  Infrasound  Monitoring 

4.1  Summary 

The  proposed  CTBT  infrasound  monitoring  network  consists  of  between  50  and  60  4- 
station  microbaragraph  arrays.  Many  of  the  infrasound  stations  will  be  co-located  or 
adjacent  to  seismic  systems  and  work  in  concert.  Each  infrasound  station  is  intended  to 
consist  of  a  broadband  microbaragraph  equipped  with  several  hundred  meters  of  noise 
reduction  hose.  The  permeable  hose  design  replaces  Daniels  microphone  pipes  for  the 
purposes  of  spatially  averaging  wind  eddy  generated  pressure  fluctuations.  Useful 
detection  thresholds  for  infrasound  stations  will  be  directly  related  to  the  effectiveness  of 
the  noise  reduction  hose  arrays.  We  present  an  analysis  of  the  differential  equations  that 
describe  the  acoustics  of  infrasound  recording  with  a  permeable  hose  as  opposed  to  the 
discrete  set  of  coupled  equations  that  have  traditionally  been  used  to  describe  a  Daniels 
pipe.  It  is  shown  that  a  hose  may  be  characterized  by  a  characteristic  time  constant,  to,  and 
a  characteristic  length,  Xo.  The  time  constant  is  related  to  permeability  of  the  hose  and  the 
characteristic  length  is  related  to  both  flow  resistance  and  permeability'  of  the  hose.  Signal- 
to-noise  improvement  is  directly  proportional  to  the  characteristic  length  of  the  hose.  The 
low  pass  filer  comer  frequency  of  the  system  is  determined  by  the  characteristic  time. 
Wavelengths  of  the  pressure  field  shorter  than  characteristic  length  are  averaged  over  the 
length  of  the  hose.  A  finite  difference  code,  Maxhose,  is  described  that  computes  response 
of  a  simple  linear  permeable  hose.  The  finite  difference  code  is  used  to  model  both 
operational  hose  designs  as  well  as  calibration  configurations.  Simulations  of  a  single  hose 
to  atmospheric  pressure  fluctuations  are  presented  for  a  white  noise  and  a  fractal  noise 
model.  A  simple  experimental  calibration  is  described  to  measure  the  characteristic  times 
and  lengths  of  permeable  hoses.  Calibration  results  are  shown  for  commercially  available 
soaker  hose.  Typical  measured  characteristic  times  are  between  10  to  20  milliseconds, 
while  characteristic  lengths  are  between  50  and  200  meters.  Of  particular  note  are  the 
effects  of  hose  degradation  during  a  typical  San  Diego  winter  as  demonstrated  by  a 
reduction  in  characteristic  length  of  the  hose  by  a  factor  of  2.  An  operational  system 
would  have  experienced  a  comparable  degradation  of  signal-to-noise  over  time.  Simple 
calibration  systems  can  be  designed  to  track  such  hose  characteristics. 

4.2  Introduction 

The  proposed  IMS  infrasound  monitoring  system  will  contain  approximately  50  to  60 
primary  infrasonic  arrays  scattered  around  the  globe.  These  arrays,  composed  of  three  or 
more  sensors,  will  often  be  co-located  or  near  an  IMS  seismic  station.  The  arrays  will  be 
designed  to  provide  arrival  time  and  azimuth  of  arrival  of  low-frequency  0.01  to  10  Hz 
sound  waves  propagating  in  the  atmosphere.  While  the  primary  mission  of  the  infrasound 
network  will  be  to  detect  nuclear  explosions  detonated  in  the  atmosphere,  these  stations 
may  also  serve  in  an  ancillary  mission  to  help  identify  large  chemical  blasts  recorded  at 
seismic  stations.  The  final  design  of  sensors  and  operational  procedures  for  detection, 
location,  and  identification  of  infrasonic  events  are  as  yet  still  in  flux.  Uncertainties  still 
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exist  in  the  areas  of  1)  excitation  as  a  function  of  source  type  (explosion,  quarry  blast, 
meteor,  volcanic  explosion)  and  size  2)  propagation  attenuation  as  function  of  distance 
and  frequency,  and  3)  optimum  sensor  design,  and  data  processing. 

Figures  40  and  41  show  projected  detection  capabilities  in  Log  10  (Yield  in  Kt)  for 
atmospheric  explosions  by  the  proposed  IMS  infrasound  network.  These  network 
simulations  show  the  proposed  network  should  be  able  to  provide  90%  assurance  of 
detection  at  2  or  more  infrasound  stations  at  the  1  Kt  level  for  most  of  the  world.  The 
90%  detection  capability  rises  above  1  Kt  in  the  South  Pacific  for  3  or  more  detections.  In 
order  to  make  projections  such  as  Figure  41,  we  require  models  for  1)  the  excitation  of 
infrasound  as  a  function  of  frequency,  yield,  and  height  of  burst,  2)  the  attenuation  of 
infrasound  as  a  function  of  frequency  and  distance,  and  3)  the  noise  and  signal  recording 
responses  of  the  recording  systems  as  a  function  of  frequency.  However,  our 
understanding  of  all  three  of  these  important  factors  are  far  from  complete.  Critical  to  the 
ability  of  the  proposed  networks  to  monitor  near  the  1  Kt  threshold  are  the  permeable 
hose  noise  reduction  systems.  Detection  thresholds  are  inversely  proportional  to  the 
signal-to-noise  improvements  expected  by  the  hose  systems. 


Figure  40.  Contours  of  detection  threshold  in  Log(Yield  in  Kt)  at  90%  probability  for 
2  infrasound  detections  at  stations  of  the  proposed  IMS.  Event  scaling, 
attenuation  relations,  and  noise  levels  based  on  Blandford  and  Clauter 
(1995)  have  been  used.  Thresholds  are  below  1  Kt  nearly  everywhere. 


36 


Figure  41 .  Contours  of  detection  threshold  in  Log  (Yield  in  Kt)  at  90%  probability  for 
3  infrasound  detections  at  stations  of  the  proposed  IMS.  Event  scaling, 
attenuation  relations,  and  noise  levels  based  on  Blandford  and  Clauter 
(1995)  have  been  used.  Thresholds  reach  2  Kt  in  the  South  Pacific. 

43  Theory  of  Noise  Reduction  Hoses  and  Their  Calibration 

The  concept  behind  the  use  of  a  noise  reduction  hose  relies  on  the  fact  that  wind 
generated  turbulent  eddies  have  shorter  apparent  wavelengths  than  an  infrasound  signal 
arriving  from  great  distance  (Cook,  1971;  Mack  and  Flinn,  1971;  McDonald  et  al.  1971). 
The  hose  “averages”  the  local  atmospheric  pressure  along  its  length  and  the  incoherent 
eddies  are  “averaged”  out. 


Figure  42.  Wind  generated  eddies  create  pressure  fluctuations  and  hence  noise.  The 
permeable  hose  acts  to  spatially  "average"  the  pressure  fluctuations. 


The  original  Daniels  pipe  (Daniels  1959)  was  conceived  as  a  discrete  set  of  inlet  ports 
along  a  length  of  pipe.  Grover  (1971)  presents  results  for  pipe  arrays  with  hypodermic 
inlet  ports,  while  Burridge  (1971)  presents  a  numerical  solution  to  this  discrete  system  of 
inlet  ports,  with  a  propagator  method  (a  system  of  coupled  equations).  However,  with  the 
advent  of  soaker  hose,  the  system  is  analogous  to  a  transmission  line  or  antennae  problem 
and  requires  the  solution  of  coupled  partial  differential  equations. 
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Figure  43.  A  hose  element:  length  dx,  radius  a,  flow  in  the  x  direction  qx,  flow  in/out  qa. 

We  consider  an  element  of  permeable  hose  of  length  dx  and  radius  a.  From 
conservation  of  mass  we  write, 

m{x,  t )  =  fix,  t)  Adx  =  q{x,  t)  =  qa  (x,  t)  -  (qx  ( x +dx,t)~  qx  (xj)), 

where  x  is  the  distance  along  the  hose,  t  is  time,  m  is  the  mass  of  air  in  a  length  of  hose 
dx,  p(x,t)=p(x,t)/(RT)  is  air  density  at  temperature  T,  A  is  the  cross  sectional  area  of  the 
hose,  q(x,t)  is  the  total  mass  flow  in  the  hose,  qa(x,t)  is  the  mass  flow  into  the  hose  from 
outside  per  unit  length,  and  qx(x,t)  is  the  lateral  mass  flow  along  the  hose  in  the  positive  x 
direction.  We  assume  that  air  diffuses  into  and  out  of  the  hose  proportional  to  the 
pressure  difference,  (p^x,t)-p(x,t»,  where  pa(x,t)  is  the  pressure  of  the  atmosphere  outside 
the  hose  and  p(x,t)  is  the  pressure  inside  the  hose, 

qx  =  pAv  =  2;zae(pa  (x,  t))  -  p(x,  t)dx , 


where  e  is  a  diffusion  constant  per  unit  length  related  to  the  permeability  of  the  hose  with 
radius,  a.  We  next  ignore  inertial  effects  and  assume  that  flow  along  the  hose  is  steady 
state  and  proportional  to  the  gradient  of  pressure  along  the  hose, 

qx  =  pAv  =  (1  /  v)dp{x,t)dx. 


where  v  is  the  hose  flow  resistance  per  unit  length  and  v  is  the  flow  velocity.  Grover 
(1971),  gives  excellent  evidence  that  flow  in  noise  reduction  pipes  can  be  described  by 
simple  Poiseuille  flow  where  the  hose  resistance  is  proportional  to  viscosity, 
v  =8tj/  pi  jt!  of  ,  where  q  is  the  viscosity  of  air.  We  combine  these  equations  to  write  a 
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partial  differential  equation  for  the  pressure  along  the  hose  in  terms  of  the  external 
atmospheric  pressure, 

2jcae  ( pa(x,t)-p(xj))-(l /  v)^^--p(x,t)A/(RT)  =  0. 

o  x 


This  equation  describes  the  hose  system  driven  by  the  atmospheric  pressure  along  its 
length.  We  transform  from  the  time-length  domain,  (x,t),  to  the  frequency  and 
wavenumber  domain,  (©,k). 


•  m 

p(x,  t )  =  p(o),  k)  exp(-/d$)  exp(-/'  Kx)dcod  k  . 

w  • 


Wre  can  write  the  solution  in  the  frequency-wavenumber  domain. 


p{(0,x)  = - - — PMr), 

(l+iW/Oc+i—f) 

Ko 


p{co,  x)  =  f - - - 7 - -pa  (<o,  k )  exp{-iKx)dK . 

J  (l  +  i(D/C)0+(~)2) 

K 0 


Now  we  observe  that  a  far-field  signal  has  the  form  pa(G),K)  -  S(x)ps(o)),  so  the  signal 
pressure  in  the  hose  is  simply 


1 

PxrnJp)  =  — — -  Psi®)  - 

^  0+16)  f&0) 

The  typical  soaker  hose  has  a  characteristic  time  constant,  to  =  2rt/©o  =  2na/{2RTe)  <  0.2 
seconds,  so  the  pressure  signal  inside  the  hose  is  a  faithful  representation  of  the  signal 
outside  the  hose  for  infrasound  frequencies  of  interest.  This  characteristic  time  is 
determined  by  the  permeability  of  the  hose.  A  very  permeable  hose  has  a  short  time 
constant  while  an  impermeable  hose  has  a  long  time  constant.  If  we  write  the  atmospheric 
noise  as  pa{co,x)  =  pR{co,K)  then  the  observed  noise  signal  is  given  by  a  Fourier  integral 
over  the  wavenumber  spectrum  of  the  atmospheric  noise, 

PnoisM  x)  =  f - - - - - PM  fc)exp(-ifcc)dK, 

J  (]+im/a>6  +  (— )2) 

Ko 

This  integral  in  general  cannot  be  evaluated  in  closed  analytic  form;  however,  by 
inspection  we  can  note  several  interesting  features.  The  incoherent  atmospheric  wavenumbers 
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are  attenuated  proportional  to  {—Y  where  ?r0  =  2mev=  47tAv! ( RT)l t0=(2tcI A0)~  is 

Ka 

the  characteristic  wavenumber  for  the  hose.  The  smaller  the  characteristic  wavenumber, 
the  more  attenuated  will  be  the  incoherent  wavenumbers  of  the  noise  field.  Note  that  the 
characteristic  wavenumber  for  the  hose  is  inversely  proportional  to  the  characteristic  time 
constant.  Therefore,  a  long  time  constant  means  more  noise  attenuation  for  a  fixed  length 
of  hose.  Typical  soaker  hoses  have  characteristic  lengths,  1c,  of  about  50-200  meters,  but 
we  have  found  that  time  constants  and  hence  characteristic  lengths  and  the  potential  noise 
reduction  change  with  time  as  the  hose  ages. 


44  L  Test  Chamber 

Pressure 


Figure  44.  A  simple  and  inexpensive  calibration  configuration. 

The  diagram  above  (Figure  44)  shows  an  experimental  setup  for  measuring  time 
constants  of  a  soaker  hose.  The  soaker  hose  is  placed  in  the  44  Liter  volume  test  chamber. 
The  test  chamber  volume  is  then  altered  with  a  12  cc  step  function  volume  decrease, 
producing  a  change  in  pressure  of  about  300  microbars.  Pressure  in  the  soaker  hose  is 
recorded  using  a  differential  pressure  gauge  (microbaragraph).  The  experiment  is  repeated 
with  different  lengths  of  hose  and  different  amplitude  step  functions  to  test  for  linearity. 
Pressure  transients  measured  with  this  setup  are  shown  in  Figure  45  for  a  virgin  soaker 
hose  and  the  same  brand  of  hose  left  outside  for  six  months  during  a  typical  San  Diego 
winter.  The  time  constant  for  tins  aged  hose  has  decreased  by  nearly  a  factor  of  two. 
Hence,  noise  reduction  will  have  decreased  by  nearly  the  same  amount. 
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5/8  Inch  Soaker  Hose  -  RiseTime  Measurements 


Figure  45.  Shows  rise  time  calibration  of  two  2.082  m  long  hoses.  The  ’’aged"  hose  had 
been  exposed  to  the  elements  for  about  6  months  during  a  typical  San  Diego 
winter.  The  rise  time  and  hence  the  characteristic  length  of  the  aged  hose  is 
about  50%  that  of  the  virgin  hose. 

Other  diagnostic  tests  are  possible  for  hose  systems.  The  diagram  below  (Figure  46) 
shows  one  such  simple  test.  The  results  are  shown  for  two  identical  lengths  (about  10  m) 
of  virgin  and  aged  soaker  hose  (Figure  47).  The  change  in  response  of  the  two  hoses  is 
clearly  evident.  In  particular,  the  aged  hose  transmits  nearly  half  of  the  pressure  amplitude 
injected  compared  to  the  virgin  hose.  This  again  demonstrates  the  importance  of 
controlling  permeability  of  the  hose  system.  This  experimental  test  setup  is  not  as  easily 
analyzed  as  the  previously  described  arrangement,  but  it  is  diagnostic.  There  is  a  tendency 
for  a  damped  Helmholtz  oscillation  to  occur.  With  shorter  lengths  of  hose  this  Helmholtz 
oscillator  can  be  used  to  measure  other  properties  of  the  hose  such  as  the  flow  resistance 
(damping). 
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Figure  46. 


L  Test  Chamber 

jw* 

A  simple,  inexpensive  diagnostic  configuration.  The  analysis  is  not  as 
simple,  but  the  test  is  sensitive  to  both  the  time  and  length  constants  of  the 
permeable  hose. 


5S  Inch  Soaker  Hose  -  Impulse  Test 


Figure  47.  Two  short  sections  of  5/8"  permeable  hose  were  tested  using  the  end-on  test 
geometry  described  in  the  text.  The  virgin  sample  shows  a  longer  time 
constant  and  a  higher  transmittance  than  the  hose  exposed  to  the  elements 
for  about  6  months.  Increased  permeability  of  the  aged  hose  has  shortened 
the  characteristic  length  and  hence  the  transmittance  of  the  hose  by  about  a 
factor  of  2  at  long  periods.  A  small  Helmholtz  oscillation  is  evident  in  the 
virgin  hose  but  was  not  excited  under  identical  condition  for  the  aged  hose. 
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4.4  Analysis  of  a  Helmholtz  Oscillator 

If  a  manifold  is  connected  to  a  short  non-penneable  hose  and  a  step  function  in 
pressure  is  applied  to  the  manifold,  it  is  possible  that  a  Helmholtz  oscillator  may  be  set  up. 
This  condition  occurs  when  the  slug  of  air  within  the  hose  protrudes  into  the  reservoir 
without  immediate  mixing  and  the  reservoir  acts  like  a  spring.  The  slug  of  air  acts  as  a 
mass  connected  to  the  spring  and  a  damped  simple  harmonic  oscillator  is  set  in  motion. 
The  viscous  flow  of  air  back  and  forth  in  the  hose  causes  damping  proportional  to  the  flow- 
velocity  and  the  damped  simple  harmonic  oscillator  can  be  modeled  to  "measure”  the 
damping  constants  used  in  the  previous  analysis. 


L 

Vr  ' 

Vi 

m  =  p  A  L 

L 

Figure  48.  Simple  diagram  of  the  principal  part  of  a  two-ended  Helmholtz  oscillator. 

The  geometry  of  the  oscillator  is  shown  above  (Figure  48).  The  air  mass  in  the  length  of 
hose  is  given  by  m  =  pAL.  The  displacement  of  the  air  mass  in  the  x  direction  is  u(t),  and 
we  write  Newton's  first  law. 


mii  =  pALii  =  -L£uJrA(Pl-Pr), 


where  Pr  and  Pi  are  the  right  and  left  manifold  pressures.  If  the  slug  of  air  protrudes  into 
the  manifold  with  displacement  u,  and  does  not  mix  with  the  air  in  the  manifold,  then  the 
effective  volume  of  the  right  and  left  manifolds  is  Vr  =  Vor  -  Au  and  Vi  =  Voi  +  Au,  where 
V,o  and  Vio  are  the  original  right  and  left  manifold  volumes.  The  right  and  left  manifold 

.  ,  SPr  ~  SVr  .  SP,  sv, 

pressures  are  given  by  — - 
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and 
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With  some  algebra,  we  have  the 
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equation  for  a  simple  harmonic  oscillator. 


pA.Lu-L^u+  A(P}  -  Pr)  =  L£u-  A2P0u 
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with  natural  frequency  Q)0  =  ^ART^l /V^-^  1/  V01) l L  —  ^jA2P0(l/  V0r-r  1 1 V01)! m  and 
damping  constant  j3  =  v  A/ 2  =  g /(2Ap) .  The  natural  period  is  proportional  to  the 
square  root  of  the  manifold  volumes  and  the  square  root  of  the  hose  length.  It  is  desirable 
to  keep  manifold  volumes  small  and  hence  keep  parasitic  Helmholtz  oscillator  frequencies 
in  the  hose  array  above  the  infrasound  frequencies  of  interest.  Likewise,  it  is  desirable  to 
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avoid  large  changes  in  the  cross  sectional  area  of  the  operational  hose-manifold  system  to 
avoid  the  Helmholtz  instability.  If  we  assume  Poiseuille  flow  in  the  hose  then  the  damping 

v'  ^  4tj 

constant  only  depends  upon  air  viscosity,  air  density,  and  hose  radius,  0  = 


The 


damping  constant  does  not  depend  upon  the  length  of  the  hose.  These  relations  can  be 
used  to  determine  damping  characteristics  and  hence  drag  terms  as  shown  in  Figure  49 
below. 


Figure  49. 


HoJmho&r  Oscflta vot  «  Notip-croncabte  Hose 


A  0.32  m  section  of  5/8”  non-permeable  hose  was  set  up  as  a  Helmholtz 
oscillator  and  subjected  to  a  0.5  millibar  step  function  change  in  pressure  of 
one  manifold.  Crosses  are  data,  the  dotted  Hne  is  a  damped  sinusoidal  fit  to 
the  data.  The  damping  constant,  x  =  l/£  =  0.37+/-0.01  sec,  is  a  measure  of 
flow  resistance  in  the  hose  and  is  consistent  with  Poiseuille  flow. 


4.5  Finite  Difference  Modeling,  Maxhose 

We  have  developed  an  implicit  finite  difference  code,  Maxhose,  to  simulate  the 
response  of  a  permeable  hose  and  manifold  system.  Our  goal  was  to  model  calibration 
configurations  such  as  those  shown  in  Figures  46  and  47  as  well  as  the  response  of  a 
system  to  an  arbitrary  atmospheric  pressure  field  specified  as  a  function  of  position  and 
time.  We  start  with  the  partial  differential  equation  for  pressure  in  the  hose  as  derived 
above. 


#*>*)  =  f^f) 


A) 


and  define  cl  =  (RT/Av)  and  c2  =  (RT/A)  2  %  a  e.  We  discretize  the  space  and  time 
variables,  xj  =  xj-i  +  dxj  and  tn  =  t„-i  +  dt*  and  define  a  pressure  vector  of  j'th  hose 
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elements,  p/n)  for  time  t*  and  pj(lrf)  for  time  W  Pressure,  pj,  is  the  pressure  of  the  volume 
between  xj  =  Xj-i  +  dx,  with  cross  sectional  area  Aj  =  x  a/  and  volume  Vj  =  Aj  dxj.  We 
follow  the  Crank-Nicolson  imphcit  method  (Smith  1978;  Davis  1986)  to  derive  a  set  of 
second  order  finite  difference  equations  for  the  pressure  updates  at  time  step  n+1  from  the 
pressures  at  time  step  n.  The  system  of  equations  that  must  be  solved  is, 

pTn  -  1 2 1  ‘&;Xd«>  -  W' + ft?) 

+(cjdtJ2)^ 

pf  +  (c‘dt„  1 2! &:)(/.<;>  -  2p^  + j£> )  -  (c’dt,  /  2)pf 
+(3'2ldtJ2)pf> 

+(c=dt„  /  2)(pi->  +  pt;>  -  2pf ) ,  forj  -  2,  N-l 

where  pj  is  the  pressure  inside  the  j'th  element  of  the  hose  and  paj  is  the  atmospheric 
pressure  outside  the  jth  element  of  the  hose.  Note  that  cl  and  c2  may  depend  upon 
position  along  the  hose  and  that  the  spatial  differences  must  be  altered  if  hose  elements  are 
not  uniform  lengths,  dxj.  We  use  second  order  approximations  for  the  spatial  and  temporal 
derivatives.  We  require  two  sets  of  boundary  conditions  at  the  ends  of  the  hose  at  xi  =  0 
and  xn  =  L.  The  first  useful  boundary  condition  is  a  no-flow  condition  where  the  spatial 
pressure  derivative  is  zero.  The  left-hand  side  no-flow  boundary  equation  becomes, 

p%!} -pT!)  =  °’  where  J=I 

and  right-hand  side  no-flow  boundary  condition 

where  j  =  N. 

The  second  set  of  useful  boundary  conditions  are  volumetric  reservoirs  or  manifolds  at 
either  end  of  a  hose.  In  the  case  of  multiple  hoses,  they  are  connected  by  manifolds  of 
finite  volume.  If  we  assume  pressure  changes  in  the  manifold  are  at  constant  temperature, 
then  we  have  the  relationship  p/(mass/V)  =  p/p  =  RT  =  po/p0  for  the  manifold.  The  change 
in  pressure  is  therefore  proportional  to  the  mass  flux, 

Sp  =  p—  =  {pf  m)q,dt  =  (RT /V)dt{U  v)  d p(x, t) !  dx. 
m 
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If  we  have  volume  reservoirs  of  volume  Vj  at  j  =  1  and/or  j  =  N,  then  using  the  Crank- 
Nicolson  scheme  we  arrive  at  the  left-hand  side  boundary  condition  equation. 
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and  the  right-hand  side  boundary  condition. 
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Equations  doatmented  up  to  this  point  assume  momentum  of  the  gas  within  the  hose  is 
insignificant  or  inertial  forces  are  not  significant.  Under  some  circumstances  this  may  not 
be  the  case,  so  we  now  derive  equations  for  the  velocities,  vj,  along  the  hose.  The  rate  of 
change  of  momentum  in  the  x-direction  in  the  hose  element  of  length,  dx,  may  be  written, 


d(mv) 

dt 


=  mv+mv~  [pAdx)v 


*  P^-A{p(x+dx)~P(x)) 

= £dxv  -  A(p(x + dx)  -  p(x)) 


= -  pvA2v  -  A[p(x + dx) — p(x)) . 


We  must  assume  that  mv«  mv  or  nonlinear  terms  (velocity  times  pressure)  will  appear  in 
the  response  of  the  hose  and  we  have. 


v=  vAv- 


'  1  ^ 


pdx  j 


{p{x+dx)~  p{x)). 


Note  that  for  constant  velocity  flow  with  a  constant  pressure  difference  along  a  section  of 

(  i  A 

A~dx  P(x))’  m<*  ^  corresponds  to  simple  Poiseuille 


hose  we  have  v  = 


V 


flow,  vpA  =  -^—(p(x+dx)~ p(x)),  and  the  flow  resistance  constants  may  be  computed 
875 dx 

from  the  air  viscosity,  £  =  8^7,  v  =  87/ k Ip! a4. 


Let  c3  =  -Av/dx  and  c4  =  -A/mass  =  -l/(pdx).  We  iet  the  discrete  velocity  vector,  vj, 
correspond  to  the  velocities  between  the  hose  elements  at  x  =  xj  at  node  j.  Therefore,  the 
velocities  and  pressures  constitute  a  staggered  grid.  Pressures  are  defined  between  xj  and 
Xj+i  while  the  velocities  are  defined  at  xj.  Again,  we  apply  the  Crank-Nicolson  formalism  to 
derive  finite  difference  equations  that  couple  velocities  and  pressures  at  time,  Vi,  and  i*. 
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dtfc3j 

“T” y 

dt,.c3. 
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dtjc4j 


Note  the  c4  coefficients  depend  upon  local  mass  and  therefore  they  are  time  dependent 
and  must  be  updated  with  each  time  step  using 


=  dtctjepp^ -pf  -p^})) 


4i 


For  left-hand  side  no-flow  boundary  conditions  we  have,v^”+;)  =  0,  forj  =  1.  For  right- 
hand  side  no-flow  boundary  conditions  we  have,  =  0,  for  N+L  If  we  have  no-flow 
boundary  conditions  on  both  ends  of  the  hose,  then  there  are  N  hose  elements  (pressures) 
and  N+l  nodes  (velocities). 

If  we  have  a  volumetric  reservoir  at  j  =  1  then  the  j  =  1  velocity  corresponds  to  the 
flow  velocity  between  the  reservoir  and  the  first  hose  element,  consistent  with  our 
previous  pressure  boundary  condition  equations,  we  have. 


If  we  have  a  volumetric  reservoir  at  j  =  N,  then  the  j  =  N-l  velocity  is  between  the  last 
(N-l)  hose  element  and  the  reservoir.  The  right-hand  side  boundary  condition  becomes. 
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dtc3. 
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For  both  right-hand  and  left-hand  side  no-flow  boundary  conditions  there  are  N  pressure 
equations  and  N+l  velocity  equations.  For  both  right-hand  side  and  left-hand  side 
volumetric  reservoirs  there  are  N  pressure  equations  and  N-l  velocity  equations.  At  each 
(n+l)  time  step,  linear  equations  are  solved,  Au(“”!>  =  Bu(n>  +  C(w<n)  -u<a))  =  b(a>,  for  the 
new  vector,  u(li+1),  of  pressures  and  velocities  W(n)  is  the  vector  of  atomospheric  pressures 
along  the  hose.  We  construct  the  pressure-velocity  vector  components;  n2j  =  pj,  U2j-i  =  Vj, 
and  W2j  =  p3j,  w?2j-i  =  0,  for  j  =  1  to  N.  C(2j-xxaj-i) =  0  and  C2j,2j  =  1  for  j  =  1,  N.  The  A  and 
B  matrices  are  very  sparse  but  not  necessarily  banded  when  multiple  coupled  hoses  are 


considered,  therefore  a  general  sparse  matrix  linear  equation  solver  is  used  rather  than  a 
solver  that  assumes  a  strictly  banded  system. 

4.6  Stability  Conditions 

While  the  Crank-Nicolson  implicit  scheme  is  generally  stable  regardless  of  time  step 
size,  we  address  one  possible  source  of  instability.  Negative  values  of  pressure  and  mass 
are  valid  solutions  to  the  equations  but  physically  meaningless.  Mass  and  pressure  are 
explicitly  positive  quantities;  therefore,  we  wish  to  avoid  large  changes  in  the  mass  of  any 
hose  or  mainfold  element,  \3m\<m  at  any  given  time  step.  The  lateral  mass  flux  across 
each  node  should  be  small  compared  to  the  mass  within  an  element  during  a  time  step, 

\dn^  =  \dt^(-^)\  <  pvol/ RT=pdx  A! RT, 


which  leads  to  a  restriction  on  the  time  step,  dt,  for  each  node. 


dt< 


p  dx  Av  Pj  (dx2AjvJ)_ _ pj_ 


R7\dp/dx\  I  ^  RT  j  jPj-Pj-ii 


dx 2 

^4  ; 


for  all  j. 


Likewise,  the  change  in  mass  of  the  right-  and  left-hand  volumetric  reservoirs  should 
remain  small  compared  to  the  mass  of  these  reservoirs  as  well  as  the  mass  of  their  adjacent 
hose  elements.  The  conditions  for  these  restrictions  on  dt  are  easily  derived.  We  write 
them  here  for  completeness. 


Wj-Ph 


and 


RT\Pj  Pj-i 


The  diffusive  flows  into  or  out  of  the  hose  elements  must  also  remain  small  for  each  time 
step. 


dm\  =  \dt2nae{pa - p)dx\  < p  vol / RT  =pdx  A/ RT, 
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for  all  j  with  e;.  >  0. 
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4.7  Hose  Calibration  Simulation 

Figure  50  shows  three  Maxhose  simulations  of  the  calibration  test  configuration 
diagramed  in  Figure  44  compared  to  virgin  hose  data  from  Figure  6.  A  proper  analysis  of 
the  calibration  data  requires  that  the  pressure  gauge  manifold  be  included  in  the  response. 
When  the  manifolds  are  included  in  the  simulation,  the  value  of  to  =  0.10  sec.  fits  the  data 
very  well. 


5 inch  Scakcr  Hose  -  Rise  Time  Measurements 


Figure  50.  Comparison  of  data  and  Maxhose  finite  difference  simulation  of  the 
calibration  test  configuration  of  Figure  44. 

4.8  Finite  Difference  Simulations  of  Wind  Noise  Response 

We  have  demonstrated  the  use  of  the  finite  difference  code.  Maxhose,  to  model  the 
hose  response  from  atmospheric  wind  noise.  An  atmospheric  pressure  history  is  specified 
as  a  function  of  time  and  position  along  the  hose,  pa(t,x),  and  the  pressure  inside  die  hose 
is  computed.  Figure  51  shows  the  results  for  pa(t,x)  specified  as  white  noise  on  a  50  m 
hose  with  characteristic  time  of  0.01  second  and  a  characteristic  length  of  100  m.  While 
this  is  not  a  realistic  wind  noise  situation,  it  illustrates  several  aspects  of  the  permeable 
hose  response.  The  numerical  experiment  was  repeated  for  several  hose  lengths  and  the 
predicted  noise  reductions  are  summarized  in  Figure  52.  RMS  noise  levels  are  reduced  by 
nearly  a  factor  40  in  the  middle  of  a  hose  that  is  longer  than  1/2  the  characteristic  length. 
Noise  reduction  saturates  at  about  1/2  the  characteristic  length.  For  such  a  hose  under 
white  noise  conditions,  noise  levels  at  the  end  of  the  hose  are  roughly  30%  higher  than 
the  levels  in  the  center  of  the  hose. 
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internal  Hose  Signal  -  White  Noise  Atmospheric  Model 
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Figure  51.  Pressure  variation  in  the  hose  as  a  function  of  time  (seconds)  and  position 
(X  in  meters)  for  random  white  noise  pressure  fluctuations  in  the 
atmosphere.  Note  that  the  hose  smoothes  out  the  pressures  along  the  hose. 


RMS  Interna!  Hose  Pressure  vs  Hose  Length 


Figure  52.  RMS  pressure  fluctuations  normalized  to  external  atmospheric  RMS 
fluctuations  at  the  end  of  the  hose  and  in  the  middle  of  the  hose  versus 
different  hose  lengths  for  a  fixed  characteristic  length  hose.  Note  that  noise 
reduction  is  best  in  the  middle  of  the  hose. 
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Under  more  realistic  wind  conditions,  atmospheric  pressure  fluctuations  exhibit 
spatial  correlation.  Given  an  arbitrary  correlation  function  or  power  spectrum  we  can 
synthesize  a  realization  of  the  atmospheric  noise,  hi  the  absence  of  a  deterministic  wind 
model,  we  have  postulated  that  pa(t,x)  appears  as  Brownian  or  fractal  noise  with  a  power 
spectrum  that  is  inversely  proportional  to  frequency  and  wavenumber,  f*a  and  k‘b,  where  a 
and  b  are  specified  constants.  Figure  53  shows  the  results  of  such  a  Monte  Carlo 
simulation  with  a  =  b  =  2.  The  same  hose  characteristics  were  used  as  in  previous 
simulations.  An  ambient  atmospheric  pressure  noise  source  is  shown  with  a  power 
spectrum  approximately  proportional  to  k'2  and  f2  is  shown  on  the  left  and  the  internal 
hose  pressure  fluctuations  are  shown  on  the  right.  Note  that  the  hose  has  spatially 
averaged  the  pressure  field.  High  frequencies  and  spatial  wavenumbers  have  been 
reduced.  However,  low  spatial  wavenumbers  are  still  evident  in  the  internal  hose  signal 
(Figure  54). 
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Figure  53.  An  atmospheric  pressure  history7  as  a  function  of  space  and  time  on  the  left 
and  the  resulting  pressure  fluctuation  signal  in  the  hose  on  the  right.  The 
50  m  hose  reduces  the  atmospheric  noise  by  averaging  it  over  the  length  of 
the  hose  with  characteristic  length  of  100  m. 
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„  Extemsf  and  Interne!  Pressure  -  Mkklle  of  Hose 


Time 


Figure  54.  Pressure  fluctuations  inside  and  outside  the  middle  of  the  hose  simulation  of 
Figure  53. 

4.9  Conclusions 

Noise  reduction  is  critical  for  the  successful  utilization  of  microbaragraph  arrays  in  a 
future  CTBT  system.  Current  plans  for  the  IMS  anticipate  signal  to  noise  increases  from 
multiple  hose  arrays  on  the  order  of  100  for  common  wind  conditions. 

We  present  an  analysis  of  the  partial  differential  equations  that  govern  the  response  of 
pressure  within  a  permeable  hose  to  external  atmospheric  pressures.  We  show  that  a  hose 
may  be  characterized  by  a  time  constant  and  characteristic  length.  The  characteristic  time 
constant  is  determined  by  the  permeability  of  the  hose.  The  characteristic  length  is 
determined  by  the  permeability  of  the  hose  and  flow  resistance  along  the  hose.  Methods 
for  measurement  of  the  time  and  length  constant  are  proposed. 

Noise  signals  in  the  hose  will  be  attenuated  for  pressure  fluctuations  in  the  atmosphere 
with  wavelengths  shorter  than  the  characteristic  hose  length.  Infrasound  signals  with 
apparent  wavelengths  longer  than  the  hose  will  not  be  so  attenuated  and  the  hose  affords  a 
net  signal-to-noise  ratio  gain.  Control  of  permeability  of  the  hose  is  critical  to  the  signal- 
to-noise  ratio  gains  that  can  be  realized.  Measurements  of  the  time  constants  of  hoses  have 
shown  that  the  permeability  can  change  by  a  factor  of  two  due  to  natural  weathering  of  a 
hose.  Methods  must  be  devised  to  protect  hoses  in  the  field  and  to  diagnose  hose 
degredation  in  operational  systems. 


Conversion  of  the  partial  differential  equations  to  a  finite  difference  program, 
Maxhose,  is  described.  Maxhose  simulates  hose  pressure  response  to  an  arbitrary 
atmospheric  pressure  defined  as  a  function  of  time  and  position  along  the  hose.  This 
second  order  implicit  Crank-Nicholson  code  runs  on  UNIX  workstations  and  Windows  95 
computers.  Results  are  shown  for  simple  calibration  configurations  and  for  two  Monte 
Carlo  realizations  of  wind  noise.  If  wind  noise  is  uncorrelated  in  space  and  time,  then  noise 
reduction  for  a  single  hose  is  predicted  to  be  about  a  factor  of  40  and  maximized  in  the 
center  of  a  hose  at  least  as  long  as  1/2  the  characteristic  length.  Noise  reduction  for  a 
single  hose  in  the  case  of  fractal  noise  is  somewhat  less  and  depends  critically  upon  the 
details  of  the  atmospheric  turbulence  spectrum.  Maxhose  can  be  extended  to  simulate 
multiple  hose  arrays  and  wind  noise  with  arbitrary  correlation  structures  in  space  and  time. 
Analysis  of  this  kind  coupled  with  an  understanding  of  wind  noise  correlation  structures 
will  prove  useful  in  understanding  the  potential  for  permeable  hose  noise  reduction. 
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